% This is the function that gets optimized by optparam. The two input arguments are the contract parameters alpha and gamma.

function output = alphagamma(y)

global beta meanp nonoilgdp oilprod pmuh pmul mean sigma 

alpha = y(1);
gamma = y(2);

expprofit1          = @(z) ((1-gamma)*z -alpha + gamma*meanp).*lognpdf(z,mean,sigma);         % This is the expected profit to the company. Oil-price - Payments to country, which are given by: alpha + gamma(p - meanp)
ExpProfit_bestcase  = quad(expprofit1,0,10000);                                               % With these parameters expected profit of no expropriation - best case that will come out of this...

if ExpProfit_bestcase < 0
   output = -ExpProfit_bestcase;
   disp(['alpha = ',num2str(alpha,4),'    gamma = ',num2str(gamma,4),'Even Bestcase makes Loss'])
end
  
%% Find pstar and capA at this combination of alpha and gamma

SEValues = alphapstar1(alpha, gamma, meanp, beta, nonoilgdp, oilprod);

p1star   = SEValues(1);
p2star   = SEValues(2);
A        = SEValues(3);
B        = SEValues(4);

%% Check whether firm's Participation constraint is met

ExpProfit  = quad(expprofit1,0,p1star) + pmuh * quad(expprofit1,p1star,p2star);      % This caluculates the expected profit to the IOC from the contract (Has to be weakly positive)

if ExpProfit < 0
    output = -ExpProfit;
    disp(['alpha = ',num2str(alpha,4),'    gamma = ',num2str(gamma,4),'     objective = ',num2str(-output,4),'    p* = ',num2str(p1star,4),'    p** = ', num2str(p2star,4),'    profits = ',num2str(ExpProfit,4)])

elseif ExpProfit >= 0
    
    %% utility of country
    
    output1 = quad(@(w)VhInt(w,p1star,p2star,A,B,alpha,gamma),0,p1star)+...
         pmuh*quad(@(w)VhInt(w,p1star,p2star,A,B,alpha,gamma),p1star,p2star)+...
         pmul*quad(@(z)Ve_l(z),p1star,10000)+...
         pmuh*quad(@(z)Ve_h(z),p2star,10000);
    
    disp(['alpha = ',num2str(alpha,6),'    gamma = ',num2str(gamma,6),'     objective = ',num2str(output1,6),'    p* = ',num2str(p1star,6),'    p** = ', num2str(p2star,6),'    profits = ',num2str(ExpProfit,6)])
    
    output  = -output1;
      
end

